Long Range Free Bosonic Models in Block Decimation Notation: Applications and 

Entanglement 
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We study the effect of long range particle exchange in bosonic arrangements. We show that by 
combining the solution of the Heisenberg equations of motion with matrix product state represen- 
tation it is possible to investigate the dynamics as well as the ground state while including particle 
exchange beyond next-neighbours sites. These ideas are then applied to study the emergence of 
entanglement as a result of scattering in boson chains. We propose a scheme to generate highly 
entangled multi-particle states that exploits collision as a powerful entangling mechanism. 



Experimental advancements in low temperature physics 
increasingly improve the control and handling of Bose 
Einstein condensates and other highly correlated states 
of matter where quantum mechanics phenomena emerge 
in their more elemental form. As a consequence, the 
gap between experiment and first principle analysis via 
clever numerical methodologies is continuously narrow- 
ing down. Powerful numerical methods as density matrix 
renormalization group and time evolving block decima- 
tion (TEBD) [l| allow to simulate systems of considerable 
complexity and size, but in most cases both methods are 
limited to short range interactions among neighbour ele- 
ments. Therefore, developing practical methods to carry 
out numerics including interactions at long scope is very 
important, as such mechanisms are integral elements of 
many physical systems, e.g., when Coulomb forces are 
involved. As it is well known, a lot of theoretical and 
experimental work has been done seeking to understand 
quantum systems and intensive research is currently car- 
ried out in many diverse fields @, 0, 0, 0, El ■ In this 
way motivated, here we explore practical alternatives for 
the study of quantum systems. We complement this dis- 
cussion by using our findings to explore the quantum 
behaviour of ID boson arrangements. We also focus on 
how much entanglement is generated after bosons collide 
in a chain. Although it has been shown that interacting 
wave packets get entangled after collision [j| [l(| , most of 
the studies on this subject focus on systems made of just 
two particles, hence applications in many-body configu- 
rations constitute an important investigation. 
Let us consider a ID arrangement of N sites in which 
boson exchange can take place at all scales. The Hamil- 
tonian is given by 
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Rk,i determines the intensity of the hopping between 
sites k and I while the creation and annihilation oper- 
ators obey the usual commuting rules [aj,ol] = 5 l k and 



[a k ,ai\ = 
(quadratic) 



= 0. Hamiltonian (JTJ is integrable 
and the Heisenberg equations of motion for 



the creation operators yield a complete set of differential 
equations of the form 
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subject to the initial condition a k (t = 0) = at. We recall 
The time dependent state ket is given 



a k = e~ itH a\e %tH . 
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nA, in such a way that the 



n k 's determine the initial configuration of bosons. Simi- 
larly, the ground state can be obtained by diagonalizing 
Hamiltonian |T]). Such ground state can be written as 
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where coefficients c k are the components of the ground 
eigenvector of matrix R, and M is the total number of 
bosons. While latter expression provides an analytical 
description of the ground ket, it does not constitute a 
practical representation of the state for large values of M 
and N. This holds especially true for computing highly 
non-local quantities which require one to deal with the 
full state. This is because using a local Fock basis to 
write the state demands exponentially growing resources. 
The purpose of this paper is to propose an alternative 
method to overcome this difficulty. Let us operate on 
the state using unitary operations with the intention of 
simplifying as much as possible Eq. ([3]). First, in order 
to make all coefficients real, we operate locally on ev- 
ery site using e~ iea ^ ak , where 9 k represents the phase of 



Cfc. 



This produces at 



Second, we operate on 



pairs of consecutive modes, {ajL i> ci k }, starting with k = 
where Q k = ■h(al +1 a k - aj^+i) 



N — 1, using e 
This induces a 



wWp 6, = ±{al +1 a k -d.{a k+li 
rotation-like transformation given by 

'7,+J. cos ( — ) - n Mil ( -f ] ] 4+i + 



c k+ ia k+1 +c k a k -> [Cfc +1 cos(^ -c^sin^) 
(c k+ \ sin {^f^j + c k cos f^f)) &k- Hence, operator a\ +1 

can always be cancelled by choosing tan (%) — 

We can carry out this cancellation consecutively until just 
the first mode is left operating on the vacuum, namely, 
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FIG. 1: Average number of bosons. In both simulations the 
hopping among sites decreases inversely with the distance, 
starting with Ri.i±i = 0.3Er for first neighbours. Er is an 
energy reference. We assume that the condensate is trapped 
in a parabolic potential with Q = 0.00046i?R Q. (a) the 
initial state corresponds to two condensates separated by a 
potential barrier. This state is obtained by calculating the 
ground state of a chain for the specified parameters plus a 
potential in the middle given by Ri^ = 1000-Er, for i = 45 — 
55. The dynamics is generated by turning the barrier off. (b) 
Starting from the ground state of the Hamiltonian, a potential 
barrier like in (a) is turned on. 



\g) = -J= (a[) |0), where any overall constant has 

been absorbed into the definition of | g) . As a result of this 
decomposition, we can obtain a representation for the 
ground state by writing \g) in block decimation notation, 
also known as matrix product states (MPS) [ll[, which 
is straightforward, and then applying the inverse unitary 
operations in reverse order using the efficient method of 
Ref. Moreover, ground state evolution, for example 
when some parameters in the original Hamiltonian are 
tuned and the whole system undergoes a quench, can also 
be studied using the same technique. In such a case the 
solution of Eq. ([2]) at any given time is inserted in equa- 
tion ^ and the reduction process takes place as above. 
From now on, we will refer to this reduction operation as 
folding. 

In order to illustrate the method, simulations in chains 
of 100 sites and 100 bosons are shown in Fig. [TJ Conden- 
sate dynamics is induced by tuning a potential barrier in 
the middle, as experimentally proposed in 0]. As a con- 
sequence, different evolution patterns emerge depending 
on how and when the potential barrier is switched on and 
off. From Fig. [T] we observe that wave trains travelling 



outwards are formed moments after a boson bulk appears 
in the centre of the chain. Then, a new bulk reassembles 
and the process starts again. The wave patterns gen- 
erated in this way are consistent with the observations 
reported in (3J, where boson trains emerge as a result of 
the rich dynamics generated in the experiment. More- 
over, the computational cost involved in the simulations 
presented above is minimal compared to TEBD compu- 
tations, where the chain must be swept many times until 
ground state convergence is achieved. 
When several summations of creation operators acting 
on the vacuum are involved, as for instance in systems 
originally arranged with bosons on different positions, 
state folding can also be implemented. For the shake of 
simplicity, let us focus on just two summations so that 
the state can be written as, 



J 



\ 



|0>, 



(4) 



/ 



where M^a account for the corresponding number of 
bosons. Si can be folded using the standard technique, 
taking into account that the coefficients in S2 are also 
affected in the process. After this, we can fold S2, but 
this time folding a\ in S2 would inevitably unfold Si. 
Therefore, the folded state acquires the following form, 



1.9) = (fit) e -«i4i (fit) |0> 



(5) 



In order to write \g) in block decimation notation, we first 
write \Mi,0, ...,0) and apply e"*^ 1 . When the state 
obtained in this way is written with explicit reference to 
the local coordinates of the first position, we can apply 

(o^ij , giving the following result [H,E|> 
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From which the canonical coefficients of state \g) can 



be directly obtained, namely, r'^" 
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and 



A'W = y/{j + Af 2 )!A 7 11 + normalization. It is impor- 
tant to note that this operation, while not unitary, does 
not modify the original configuration of Schmidt vectors 
across the chain. To see this, consider first the Schmidt 
vectors corresponding to the first site, {^i)' 1 '}. Because 
Hamiltonian ([1]) preserves the total amount of particles, 
these vectors each have a definite number of bosons as- 
sociated with them. Consequently, applying d} lifts the 
state occupation by one. These lifted kcts are orthogonal 
and therefore valid Schmidt vectors of the whole updated 
state. On the other hand, any other Schmidt decomposi- 
tion of the system is made of only one vector to the left 
and one vector to the right, since there is no boson at all 
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between sites 3 and N. As a result, lifting the local basis 
in the first site cannot alter the topology of the original 
decomposition. 

In order to illustrate how these ideas can be applied to 
solve challenging problems, we first introduce our dynam- 
ical model analytically so as to synthesize some useful 
results. Consider the case where M bosons are sent from 
one end of the chain to the other [2|, |9( under the pres- 
ence of a potential barrier in the central part. Evolution 
is therefore given by df f |0). We assume that the dy- 
namics is described by a matrix R = J x + ee~P J * , where 
Jx,y,z & r e the standard angular momentum operators. In 
this way, J x induces transmission of bosons from one end 
of the chain to the other while the term proportional to 
e represents a perturbation symmetrically localized over 
the central part of the chain. R acts on the set of eigen- 
vectors of J z , which are associated to the original op- 
erators according to \j, m) — > dj_ m+1 , where j — 
Similarly, the evolution of operator &\ can be studied by 
following the time evolution of ket \j, j): 

m^-K))^e-^+^^\j,j) (7) 

« \j,-j)-ie r dte-^- t)J *e-P p >e- uJ *\j,j), 
Jo 

In latter equation we have expanded the time evolution 
operator to first order in e. The integral can also be 
written as 

f 7 dte- l{7T - t] - h e- (i3 * ''e l{7T ~ t)S * e~ l7rJ *\j,j) 
Jo 

= f dte- picos< - t)J *- sMt)J « )2 \j,-j) (8) 
Jo 

V4^5 J- x Jo 

In the last line we have used the well known identity 
f°° dxe~^ ax ' +bx ^ — */%e%z . We can calculate the ac- 

J — oo V a 

tion of angular momentum operators in the integral 
above using Schwinger's model of angular momentum. 
To this end, we define two sets of creation and annihi- 
lation operators d + ,d^_ and d_,d^, in terms of which 
Jz — ^(<^a + — a_aJ), J + = a^_a_, J_ = a^_a + and 

\j,m) = — % = 10). In this way integral JHJ) be- 

y/(j+m)\(j-m)\ 

comes, 
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FIG. 2: Entanglement between the ends vs. the normal- 
ized intensity of the potential barrier, fi = ee - ^ (arbitrary 
energy units), in a chain of 100 sites. The right axis coordi- 
nates correspond to the solid line, the amount of bosons on 
the ends divided by M. In this case curves corresponding 
to different number of bosons are very similar and therefore 
they all collapse in one single line. Inset. Intensity of en- 
tanglement maxima in the main figure as a function of the 
number of bosons. 

power series with little loss in accuracy when the negative 
exponential in the first integral suppresses the contribu- 
tion of the second integral for relatively high values of 
y. This can be guaranteed as long as f3 <C |, which ac- 
counts for a perturbation operating on a section of the 
chain much smaller than the system size. With this ap- 
proximation equation ([9]) gives, 




r (\^ + ^) J ^^vMJ-k), (10) 



which can be integrated using the asymptotic properties 
of the Bessel functions. So, replacing state vectors by 
operators we obtain, 





where we have performed the variable change x = jy. 
We then use the binomial expansion and obtain a sum of 
integrals. We can expand the trigonometric functions in 



Consequently, for small /3, operator di will generate a 
particle configuration exponentially localized around the 
end opposite to the position where particles were initially 
localized, that is, particles are efficiently transferred from 
one end to the other. On the other hand, for larger (3 the 
main contribution will come from a\ as well as from ajy. 
Hence particles will occupy both ends with little spread- 
ing over intermediate sites. Physically, large (3 means the 
perturbation remains tightly localized in the centre of the 
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chain. In this case one can think that when the parti- 
cles get across the centre they have a well defined mo- 
mentum and the perturbation simply acts as a thin wall 
that causes reflection without altering the wave packet 
shape. As a result, particles reflected preserve the co- 
herence necessary to be drifted back into their original 
chain terminal. This phenomenon can be used to effi- 
ciently generate entanglement among colliding particles. 
Ideally boson packets initially prepared far away from 
each other in a separable state interact in the middle of 
the chain via a local potential proportional to the num- 
ber of boson on the central positions. Some time after 
the interaction, the majority of bosons turn up on the 
original positions, but this time these bosons display long 
range entanglement. This process has been simulated us- 
ing the two-sum-folding technique previously introduced. 
So, we solve equation {2j and insert this solution in the 
expression that determine the initial state as previously 
specified. Then, for a time equal to half a period we set 
the quantum state in MPS representation and find the 
entanglement between the ends of the chain as measured 
by log- negativity [!, @j- Here we study chains with N 
even and therefore a potential barrier is considered in two 
central sites, but in chains with N odd a highly localized 
perturbation can be modelled taking one single central 
site. From Fig. [5] it can be seen that particle collec- 
tion is highly efficient, with more than 99% of particles 
collected on the terminals after the interaction. Simi- 
larly, entanglement generation is maximally efficient in 
the range 0.5 < n/N < 1.5 (see Figf2]caption for the def- 
inition of fi) while the maximum amount of entanglement 
grows logarithmically as a function of the total number 
of bosons. This logarithmic behaviour can be studied 
applying folding techniques. To this end, consider the 

M M 

dynamics generated by &■£ |0). Because bosons get 
entangled mainly on the central part of the chain where 
the particle packets interact, we focus on the form of a- 
operators for a time equal to a quarter period, when the 



boson packets become symmetrically spread along the 
chain and therefore di(-j) = o-n^) for a chain without 
wall. Consequently, the particle distribution adopts the 
form of Eq. ([3]), but this time the the c^'s are determined 
by the coupling constants and the strength of the barrier, 
both considered general up to symmetry. As entangle- 
ment between both halves of the chain is independent 
of unitary transformations on either half-chain block, we 
can operate on the state using the same transformations 
employed previously to fold the quantum state as long 
as we do not operate on pairs of creation operators be- 
longing to different half-chain blocks. In this case, the 
result of the reduction is trivial and can be written as 
(a\ + a2) M |0)- Entanglement can be calculated expand- 
ing this expression and noticing that as a result we ob- 
tain a genuine Schmidt decomposition with coefficients 

Ck = y 'jM -k)\k\ ■ Hence, log-negativity can be worked 
out by replacing the binomial distribution by a normal 
Gaussian. This gives -Eat ~ ^0(72(2^) + \logi (M). More- 
over, assuming that the interaction among bosons takes 
place in a time scale much smaller than the transfer time, 
we can take this result as an estimation of entanglement 
when the particles reach the chain ends. The discrep- 
ancy with the fitting constants in Fig. [3] is due to to the 
approximations involved in getting En- However, log- 
arithmic growing as well as highly efficient particle col- 
lection are positively verified by the numerical analysis. 
Finally, we would like to point out that the methodology 
presented in this work can be applied to analogous situa- 
tions in spin chains or fermion arrangements. Non linear 
effects, while less amenable to exact analytical solution, 
could be considered by solving the Gross-Pitaevskii equa- 
tion and then treating the coefficients describing the wave 
function as the Cfc's of Eq. ([3]). Additionally, the method 
can be used in combination with TEBD, for example in 
quench-like problems where repulsion is suddenly turned 
on. JR acknowledges an EPSRC-DHPA scholarship. 
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